function y = func_Km_R_nonmod(zeta, m, CONSTS)

    debug_output = false;
    plot_data = false;
    
    k0 = CONSTS.k0;
    a  = CONSTS.a;
    d  = CONSTS.d;

    y = zeros(size(zeta));
    
    for i = 1:size(zeta, 1)
        zeta_i = zeta(i);

        period = pi/(k0*a);
        spans_per_period = 18;
        order = 3; %3; %3 %1
        if (size(zeta, 1) >= 2)
            adding = 0.9*k0*zeta(2);
        else
            adding = 0.05*d;
        end
        
        q_high = 0.5*log(10)*order/(k0*(zeta_i) + adding);

        span_size = period/spans_per_period;
        
        integ1 = quadgk(@(q)(real(func_ns_nonmod(q, zeta_i, m, CONSTS))), ...
            0, 1); %1.001
        integ2 = quadgk(@(q)(imag(func_ns_nonmod(q, zeta_i, m, CONSTS))), ...
            1, inf); %1.001

%         integ1 = quad_gauss(@(q)(real(func_ns_nonmod(q, zeta_i, m, CONSTS))), ...
%             [0, 1], span_size, 6); %1.001
%         integ2 = quad_gauss(@(q)(imag(func_ns_nonmod(q, zeta_i, m, CONSTS))), ...
%             [1, q_high], span_size, 6); %1.001
        integ = integ1+1i*integ2;

        if (debug_output)
            fprintf('i=%d    %e + %ei\n', i, ...
                real(integ), imag(integ));
        end
        
        y(i) = integ;
    end
    
    if(plot_data)
        figure; plot(zeta./d, real(y), 'b.-', zeta./d, imag(y), 'r.-');
        title('K_m^{(R)}(\zeta)'); legend('Re(K_{m, reg-nonmod})', 'Im(K_{m, reg-nonmod})');
    end
    
    if(false)
        y1 = kernel_m_singular_nonmod_1_analytical(zeta, m, CONSTS);
        y2 = kernel_m_singular_nonmod_2_analytical(zeta, m, CONSTS);
        K_s = y1+y2;
        figure; plot(zeta./d, real(K_s), 'b-', zeta./d, imag(K_s), 'r-');
        title('K_m^{(S)}(\zeta)'); legend('Re(K_{m, s-nonmod})', 'Im(K_{m, s-nonmod})');
    end
   
end

